rm(list=ls())
library(tidyverse)
setwd("~/Downloads/dataverse_files(3)")

getmode <- function(v, na.omit = TRUE){
  if (na.omit == TRUE) {
    v <- na.omit(v)  
  }
  uv <- unique(v)
  tab <- tabulate(match(v, uv))
  out <- uv[tab == max(tab)]
  if(length(out) > 1){
    out <- sample(out,1)
  }
  return(out)
}

# add custom functions to extract estimates (tidy) and goodness-of-fit (glance) information
tidy.feglm <- function(x, ...) {
  s <- summary(x,  type = "sandwich")
  ret <- data.frame(
    term      = row.names(s$cm),
    estimate  = s$cm[, 1],
    std.error  = s$cm[, 2]
  )
  ret
}

glance.feglm <- function(x, ...) {
  ret <- data.frame(
    Null.Dev = x$null.deviance,
    Deviance = x$deviance,
    N   = x$nobs[4],
    `Committee FE` = as.integer(length(x$nms.fe[[1]])),
    `Congress FE` = as.integer(length(x$nms.fe[[2]])))
  ret
}


matched_dat<- read_csv("Data/ParkGrandstandingMatchedData.csv")

cor.test(matched_dat$num_citations, matched_dat$avg_grandstand)

library(fixest)
num_cite_grand <- fenegbin(num_citations ~ party_control*avg_grandstand | Committee_short, data = matched_dat)

matched_dat$cites_science <- 0
matched_dat$cites_science[matched_dat$num_citations > 0 ] <- 1
logit_cite_grand <- feglm(cites_science ~ party_control*avg_grandstand | Committee_short, data = matched_dat, family = binomial(link = "logit"))

library(marginaleffects)

negbin_plot <- plot_predictions(num_cite_grand, condition= c("avg_grandstand" ,"party_control")) + 
  theme_minimal() + 
  scale_color_manual("", values=c("#9A3E25", "#156B90")) + 
  scale_fill_manual("", values=c("#9A3E25", "#156B90")) + 
  ylab("Predicted citations to science") + 
  xlab("Hearing Average Grandstanding Score (Park 2021)") +
  theme(legend.position = "bottom")

logit_plot <- plot_predictions(logit_cite_grand, condition= c("avg_grandstand" ,"party_control")) + 
  theme_minimal() + 
  scale_color_manual("", values=c("#9A3E25", "#156B90")) + 
  scale_fill_manual("", values=c("#9A3E25", "#156B90")) + 
  ylab("Predicted probability of citing science") + 
  xlab("Hearing Average Grandstanding Score (Park 2021)") +
  theme(legend.position = "bottom")


library(cowplot)

grand_plot <- cowplot::plot_grid(negbin_plot, logit_plot, labels = c("A", "B"))

ggsave("Output/FigS30.pdf", grand_plot, width = 9, height = 5)
